Load in Libs

library(broom)
## Warning: package 'broom' was built under R version 3.6.3
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 3.6.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.6.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages ---------------------------------------------------------------------------- tidyverse 1.3.0 --
## v tibble  3.0.3     v dplyr   1.0.0
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## v purrr   0.3.4
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'readr' was built under R version 3.6.3
## Warning: package 'dplyr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(mgcv)
## Warning: package 'mgcv' was built under R version 3.6.3
## Loading required package: nlme
## Warning: package 'nlme' was built under R version 3.6.3
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(plotly)
## Warning: package 'plotly' was built under R version 3.6.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Read in raw data from RDS.

raw_data <- readRDS("./n1_n2_cleaned_cases.rds")

final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
  rename(Facility = wrf) %>%
  mutate(Facility = recode(Facility, 
                           "NO" = "WRF A",
                           "MI" = "WRF B",
                           "CC" = "WRF C"))

Data split for main figure

#split data to layer for main plotly figure
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>% 
  select(c(date, cases_cum_clarke, new_cases_clarke, X10_day_ave_clarke, cases_per_100000_clarke)) %>%
  group_by(date) %>% summarise_if(is.numeric, mean)



#specify fun colors
background_color <- "#7570B3"
ten_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')
#creates the two panels of the main figure
        p1 <- only_background %>%
              plotly::plot_ly() %>%
              plotly::add_trace(x = ~date, y = ~new_cases_clarke, 
                                type = "bar", 
                                hoverinfo = "text",
                                text = ~paste('</br> Date: ', date,
                                                     '</br> Daily Cases: ', new_cases_clarke),
                                alpha = 0.5,
                                name = "Daily Reported Cases",
                                color = background_color,
                                colors = background_color,
                                showlegend = FALSE) %>%
            layout(yaxis = list(title = "Athens Daily Cases", range = c(0,80), showline=TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #renders the main plot layer two as ten day moving average
        p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X10_day_ave_clarke, 
                             type = "scatter",
                             mode = "lines",
                             hoverinfo = "text",
                            text = ~paste('</br> Date: ', date,
                                                     '</br> Ten-Day Moving Average: ', X10_day_ave_clarke),
                             name = "Ten Day Moving Average Athens",
                             line = list(color = ten_day_ave_color),
                             showlegend = FALSE)
        
        #renders the main plot layer three as positive target hits
        p2 <- plotly::plot_ly() %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n1,
                                       symbol = ~Facility,
                                       marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n2,
                                       symbol = ~Facility,
                                       marker = list(color = '#D95F02', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
            layout(yaxis = list(title = "SARS CoV-2 Copies/L", 
                                 range = c(3, 8), showline = TRUE,
                                 type = "log",
                                 automargin = TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #adds LOD line
        p2 <- p2 %>% plotly::add_segments(x = as.Date("2020-03-14"), 
                                          xend = ~max(date + 10), 
                                          y = 3571.429, yend = 3571.429,
                                          opacity = 0.35,
                                          line = list(color = "black", dash = "dash")) %>%
          layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y", 
                                    text = "Limit of Detection", showarrow = FALSE))

        

        p1
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: Ignoring 1 observations
        p2
## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
#create smoothing data frames 
#n1
smooth_n1 <- only_n1 %>% select(-c(Facility)) %>% 
  group_by(date, cases_cum_clarke, new_cases_clarke, X10_day_ave_clarke, cases_per_100000_clarke) %>%
  summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
  ungroup() %>%
  mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
  mutate(target = "N1")
## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X10_day_ave_clarke' (override with `.groups` argument)
#n2
smooth_n2 <- only_n2 %>% select(-c(Facility)) %>% 
  group_by(date, cases_cum_clarke, new_cases_clarke, X10_day_ave_clarke, cases_per_100000_clarke) %>%
  summarize(sum_copy_num_L = sum(mean_total_copies)) %>%
  ungroup() %>%
  mutate(log_sum_copies_L = log10(sum_copy_num_L)) %>%
  mutate(target = "N2")
## `summarise()` regrouping output by 'date', 'cases_cum_clarke', 'new_cases_clarke', 'X10_day_ave_clarke' (override with `.groups` argument)
#loess at 0.80 span (same as ggplot default)
sumfit_1 <- loess(log_sum_copies_L ~ new_cases_clarke, data = smooth_n1, span = 0.8)
sumfit_2 <- loess(log_sum_copies_L ~ new_cases_clarke, data = smooth_n2, span = 0.8)
#generalized additive model 
sumfit_3 <- gam(smooth_n1$log_sum_copies_L ~ smooth_n1$new_cases_clarke)
sumfit_4 <- gam(smooth_n2$log_sum_copies_L ~ smooth_n2$new_cases_clarke)
#generalized linear model
sumfit_5 <- glm(log_sum_copies_L ~ new_cases_clarke, data = smooth_n1)
sumfit_6 <- glm(log_sum_copies_L ~ new_cases_clarke, data = smooth_n2)

p3 <- plotly::plot_ly() %>%
  plotly::add_trace(x = ~date, y = ~log_sum_copies_L,
                    type = "scatter",
                    mode = "markers",
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date,
                                  '</br> Copies/L: ', round(log_sum_copies_L, digits = 2)),
                    data = smooth_n1,
                    marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
  plotly::add_trace(x = ~date, y = ~log_sum_copies_L,
                    type = "scatter",
                    mode = "markers",
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date,
                                  '</br> Copies/L: ', round(log_sum_copies_L, digits = 2)),
                    data = smooth_n2,
                    marker = list(color = '#D95F02', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
  plotly::add_lines(x = ~date, y = predict(sumfit_1),
                    data = smooth_n1,
                    hoverinfo = "text",
                    text = NULL,
                    showlegend = FALSE,
                    line = list(color = '#1B9E77')) %>%
  plotly::add_lines(x = ~date, y = predict(sumfit_2),
                    data = smooth_n2,
                    hoverinfo = "text",
                    text = NULL,
                    showlegend = FALSE,
                    line = list(color = '#D95F02'))%>%
              layout(yaxis = list(title = "Log SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 automargin = TRUE))
  


p3
test12 <-
    plotly::subplot(p3,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
test12
#make a dataframe to average both genes, remove target column, group by date, then average
both_clean <- full_join(smooth_n1, smooth_n2) %>%
  select(-c(target)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup()
## Joining, by = c("date", "cases_cum_clarke", "new_cases_clarke", "X10_day_ave_clarke", "cases_per_100000_clarke", "sum_copy_num_L", "log_sum_copies_L", "target")
combofit_1 <- loess(log_sum_copies_L ~ new_cases_clarke, data = both_clean, span = 0.8)
combofit_2 <- lm(log_sum_copies_L ~ new_cases_clarke, data = both_clean)
p4 <- plotly::plot_ly() %>%
  plotly::add_trace(x = ~date, y = ~log_sum_copies_L,
                    type = "scatter",
                    mode = "markers",
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date,
                                  '</br> Copies/L: ', round(log_sum_copies_L, digits = 2)),
                    data = both_clean,
                    marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
    plotly::add_lines(x = ~date, y = predict(combofit_1),
                    data = both_clean,
                    hoverinfo = "text",
                    text = NULL,
                    showlegend = FALSE,
                    line = list(color = '#1B9E77'))%>%
              layout(yaxis = list(title = "Log SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 automargin = TRUE))
  

p4
#make a plot that stacks the averge of both genes
test1234 <-
    plotly::subplot(p4,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
test1234
#add trendlines 
#extract data from geom_smooth
#n1 extract
cc <- ggplot(smooth_n1, aes(x = date, y = log_sum_copies_L)) + 
  stat_smooth(aes(outfit=fit_n1<<-..y..), color = '#1B9E77', span = 0.8)
## Warning: Ignoring unknown aesthetics: outfit
#n2 extract
oo <- ggplot(smooth_n2, aes(x = date, y = log_sum_copies_L)) + 
  stat_smooth(aes(outfit=fit_n2<<-..y..), color = '#1B9E77', span = 0.8)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#n1
cc
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

fit_n1
##  [1]  9.599321  9.795726  9.985385 10.168217 10.344138 10.513066 10.674920
##  [8] 10.829616 10.977062 11.116207 11.246872 11.369763 11.485585 11.595042
## [15] 11.698840 11.797682 11.892275 11.983020 12.069420 12.150825 12.226582
## [22] 12.296041 12.358548 12.413454 12.460118 12.498344 12.528405 12.550586
## [29] 12.565168 12.572437 12.572675 12.566165 12.551775 12.518279 12.467361
## [36] 12.404187 12.333921 12.261729 12.192777 12.132229 12.080847 12.017299
## [43] 11.942145 11.860341 11.776843 11.696605 11.624585 11.565736 11.516511
## [50] 11.448078 11.364294 11.272870 11.181518 11.097947 11.029870 10.984996
## [57] 10.963488 10.947805 10.937376 10.932909 10.935114 10.944700 10.962375
## [64] 10.988850 11.023697 11.065042 11.113084 11.168088 11.230321 11.300047
## [71] 11.377533 11.463044 11.556458 11.657270 11.765407 11.880798 12.003373
## [78] 12.133062 12.269794 12.413498
#n2
oo
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

fit_n2
##  [1]  9.501965  9.703836  9.900223 10.091067 10.276310 10.455894 10.629761
##  [8] 10.797852 10.960110 11.116574 11.267296 11.412276 11.551515 11.685015
## [15] 11.812776 11.934799 12.051247 12.166930 12.282540 12.395966 12.505093
## [22] 12.607810 12.702005 12.785564 12.856663 12.918351 12.972096 13.017976
## [29] 13.056068 13.086451 13.109204 13.124403 13.130773 13.117455 13.086028
## [36] 13.041404 12.988497 12.932220 12.877486 12.829210 12.787393 12.727972
## [43] 12.652198 12.566285 12.476449 12.388903 12.309863 12.245544 12.194276
## [50] 12.129610 12.053696 11.972170 11.890668 11.814826 11.750281 11.702670
## [57] 11.672816 11.649440 11.631228 11.617657 11.608204 11.602344 11.599556
## [64] 11.599317 11.601726 11.607889 11.617951 11.632020 11.650207 11.672620
## [71] 11.699370 11.730564 11.766164 11.805962 11.849903 11.897934 11.950000
## [78] 12.006050 12.066028 12.129882
#assign fits to a vector
n1_trend <- fit_n1
n2_trend <- fit_n2

#extract y min and max for each
bb<- ggplot_build(cc)$data
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
bb<- as.data.frame(bb)
n1_ymin <- bb$ymin
n1_ymax <- bb$ymax

qq <- ggplot_build(oo)$data
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
qq <- as.data.frame(qq)
n2_ymin <- qq$ymin
n2_ymax <- qq$ymax

#reassign dataframes (just to be safe)
work_n1 <- smooth_n1
work_n2 <- smooth_n2

#fill in missing dates to smooth fits
work_n1 <- work_n1 %>% complete(date = seq(min(date), max(date) +9, by = "1 day"))
date_vec_n1 <- work_n1$date
work_n2 <- work_n2 %>% complete(date = seq(min(date), max(date) +9, by = "1 day"))
date_vec_n2 <- work_n2$date

#create a new smooth dataframe to layer
smooth_frame_n1 <- data.frame(date_vec_n1, n1_trend, n1_ymin, n1_ymax)
smooth_frame_n2 <- data.frame(date_vec_n2, n2_trend, n2_ymin, n2_ymax)

#plot smooth frames
p99 <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_n1, y = ~n1_trend,
                    data = smooth_frame_n1,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n1,
                                  '</br> Log Copies/L: ', round(n1_trend, digits = 2),
                                  '</br> Target: N1'),
                    line = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
plotly::add_lines(x = ~date_vec_n2, y = ~n2_trend,
                  data = smooth_frame_n2,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n2,
                                  '</br> Log Copies/L: ', round(n2_trend, digits = 2),
                                  '</br> Target: N2'),
                    line = list(color = '#D95F02', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
plotly::add_ribbons(x ~date_vec_n1, ymin = ~n1_ymin, ymax = ~n1_ymax,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n1, #leaving in case we want to change
                                  '</br> Max Log Copies/L: ', round(n1_ymax, digits = 2),
                                  '</br> Min Log Copies/L: ', round(n1_ymin, digits = 2),
                                  '</br> Target: N1'),
                    name = "",
                    line = list(color = '#1B9E77')) %>%
plotly::add_ribbons(x ~date_vec_n2, ymin = ~n2_ymin, ymax = ~n2_ymax,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_n2, #leaving in case we want to change
                                  '</br> Max Log Copies/L: ', round(n2_ymax, digits = 2),
                                  '</br> Min Log Copies/L: ', round(n2_ymin, digits = 2),
                                  '</br> Target: N2'),
                    name = "",
                    line = list(color = '#D95F02')) %>%
                layout(yaxis = list(title = "Log SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date"))

p99
smooth_extract <-
    plotly::subplot(p99,p1, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )
## Warning: Ignoring 1 observations
smooth_extract
save(smooth_extract, file = "./smooth_extract.rda")

Static options

mindate <- min(raw_data$date)
maxdate <- max(raw_data$date)
#static trend plot using loess span 0.8
trend_plot <- ggplot(smooth_n1, aes(x = date, y = log_sum_copies_L)) + 
  geom_smooth(color = '#1B9E77', span = 0.8) +
  geom_point(color = '#1B9E77') +
  geom_smooth(data = smooth_n2, color = '#D95F02', span = 0.8) +
  geom_point(data = smooth_n2, color = '#D95F02', span = 0.8)
## Warning: Ignoring unknown parameters: span
trend_plot <- trend_plot + theme_classic() + xlab("Date") + ylab("Log Copies/L") + 
  scale_x_date(limits = c(mindate, maxdate))

trend_plot 
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#static epi curve for facet wrapping
epi_curve <- ggplot(only_background, aes(x = date, y = new_cases_clarke)) + 
  geom_bar(stat = "identity", fill = "#7570B3", alpha = 0.5) +
  geom_line(aes(x = date, y = X10_day_ave_clarke), color = "#E6AB02") +
  xlab("Date") + ylab("Athens Daily Cases") +
  theme_classic()

epi_curve
## Warning: Removed 1 rows containing missing values (position_stack).
## Warning: Removed 10 row(s) containing missing values (geom_path).

smooth_stack <- gridExtra::grid.arrange(trend_plot, epi_curve, nrow =2)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 10 row(s) containing missing values (geom_path).

smooth_stack
## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
ggsave(smooth_stack, file = "./smooth_stack.png")
## Saving 7 x 5 in image
combogg <- ggplot(both_clean, aes(date, log_sum_copies_L)) + 
         geom_smooth(method = "loess", span = 0.8, color = '#1B9E77', alpha = 0.5, size = 1) +
         geom_point(size=1) +
  scale_x_date(limits = c(mindate, maxdate))

combo_gg_plotly <- ggplotly(combogg)
## `geom_smooth()` using formula 'y ~ x'
combo_gg_plotly
epi_gg_plotly <- ggplotly(epi_curve)
## Warning: Removed 1 rows containing missing values (position_stack).
combo_combo <-
    plotly::subplot(combo_gg_plotly,epi_gg_plotly, # plots to combine, top to bottom
      nrows = 2,
      heights = c(.6,.4),  # relative heights of the two plots
      shareX = TRUE,  # plots will share an X axis
      titleY = TRUE
    ) %>%
    # create a vertical "spike line" to compare data across 2 plots
    plotly::layout(
      xaxis = list(
        spikethickness = 1,
        spikedash = "dot",
        spikecolor = "black",
        spikemode = "across+marker",
        spikesnap = "cursor"
      ),
      yaxis = list(spikethickness = 0)
    )

combo_combo